Associations Between Undertrition in Children With the Standardised Precipitation Evapotranspiration Index (SPI) of municipalities in Brazil (2008-2019)

Author

Arthur R. Magalhães & Daniel Vartanian

Published

2025-01-21

Overview

This document explore possible associations between undertrition in children with the Standardised Precipitation Evapotranspiration Index (SPEI) of municipalities in Brazil durting 2008-2019 period. It is part of the Sustentarea Research and Extension Group’s project titled Global Syndemic: The Impact of Anthropogenic Climate Change on the Health and Nutrition of Children Under Five Years Old Served by Brazil’s Public Health System (SUS).

Please note that this is a preliminary report intended to support decision-making and may not include all the details of the analysis.

Methods

Approach and procedure method

The analysis will rely solely on a visual inspection of the data. If necessary, additional analysis and updates will be incorporated into this report later.

Source of data/information

The data used in this analysis have as sources:

  • The Brazilian Institute of Geography and Statistics (IBGE) Automatic Retrieval System (SIDRA), for data on GDP per capita.
  • Brazil’s Food and Nutrition Surveillance System (SISVAN) (link to SISVAN, for data on malnutrition.
  • WorldClim, for data on bioclimatic variables, which allowed us to calculate the Standardised Precipitation Evapotranspiration Index (SPEI) for municipalities in Brazil.

Data wrangling

The data wrangling and analysis followed the data science framework outlined by Hadley Wickham et al. (Wickham et al., 2023), as illustrated in Figure 1.

Figure 1: Model of the data science process created by Wickham, Çetinkaya-Runde, and Grolemund (2023).

Source: Reproduced from Wickham et al. (2023).

All the analyses are 100% reproducible and can be run again at any time. See the README file on the code repository to learn how to run them.

Setting the Enviroment

Code

Loading functions

Code
source(here::here("R", "coef_mean.R"))
source(here::here("R", "model_gamm_misfs.R"))
source(here::here("R", "plot_gamm.R"))
source(here::here("R", "utils.R"))
source(here::here("R", "utils-plots.R"))

Importing and Tidying the Data

GDP Data

Code
gdp_data <- 
  here::here("data", "PIB dos Municípios - base de dados 2010-2021.xlsx") |>
  readxl::read_excel(col_types = "text") |>
  janitor::clean_names() |>
  dplyr::select(
    ano, 
    codigo_do_municipio, 
    produto_interno_bruto_per_capita_a_precos_correntes_r_1_00
  ) |>
  dplyr::rename(
    year = ano,
    municipality_code = codigo_do_municipio,
    gdp_per_capita = 
      produto_interno_bruto_per_capita_a_precos_correntes_r_1_00,
  ) |>
  dplyr::mutate(
    year = as.integer(year),
    municipality_code = 
      municipality_code |>
      stringr::str_sub(end = -2) |>
      as.integer(),
    gdp_per_capita = as.numeric(gdp_per_capita)
  )

gdp_data

Nutrition Data

Code
nutrition_data <-
  here::here("data", "Banco_dados_malnutritio_clima - Adaptado.csv") |>
  readr::read_csv(col_types = readr::cols(.default = "c")) |>
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::everything(),
      .fns = ~ 
        iconv(.x, from = "UTF-8", to = "iso-8859-1") |>
        readr::parse_guess()
    )
  ) |>
  janitor::clean_names() |>
  dplyr::rename(
    year = ano,
    municipality_code = code_muni
  ) |>
  dplyr::mutate(
    year = as.integer(year),
    municipality_code = as.integer(municipality_code),
    cobrs = as.numeric(cobrs)
  )
#> New names:
#> • `` -> `...1`

nutrition_data

SPEI Data

Code
spei_data <- 
  here::here("data", "spei_Extreme_drought_event_municipality_year2.csv") |>
  readr::read_csv(col_types = readr::cols(.default = "c")) |>
  janitor::clean_names() |>
  dplyr::rename(
    municipality_code = code_muni,
    municipality_name = name_muni
  ) |>
  dplyr::select(
    municipality_code, 
    dplyr::starts_with("spei_12m_20")
  ) |>
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::starts_with("spei_12m_20"),
      .fns = as.numeric
    )
  ) |>
  tidyr::pivot_longer(
    cols = starts_with("spei_12m_20"),
    names_to = "year",
    values_to = "spei_12m"
  ) |>
  dplyr::mutate(
    municipality_code = 
      municipality_code |>
      substr(1, nchar(as.character(municipality_code)) - 1) |>
      as.integer(),
    year = 
      gsub("spei_12m_", "", year) |>
      as.integer(),
  )
  
spei_data

Analysis Data

Code
data <- 
  gdp_data |>
  dplyr::full_join(
    y = nutrition_data,
    by = c("year", "municipality_code")
  ) |>
  dplyr::left_join(
    spei_data,
    by = c("year", "municipality_code")
  )

data

Checking for Muncipalities With Minimum 5% Cover from SISVAN

data |>
  dplyr::pull(cobrs) %>%
  sum(. < 0.05, na.rm = TRUE)
#> [1] 29680.3019
data |>
  dplyr::pull(cobrs) %>%
  sum(. > 1, na.rm = TRUE)
#> [1] 30664.3019
data <- 
  data |>
  dplyr::filter(cobrs >= 0.05)

data
min(data$spei_12m)
#> [1] -1.751599753
max(data$spei_12m)
#> [1] 1.026890197

Generalized Linear Mixed Model (GLMM)

MBEPR – Very Short Stature for Age (Muito Baixa Estatura Para a Idade)

See https://cran.r-project.org/web/packages/glmmTMB/vignettes/troubleshooting.html to learn more about the model warnings.

mbepr_model <- glmmTMB::glmmTMB(
  mbepr ~ spei_12m + gdp_per_capita + (1 | year),
  data = data,
  family = glmmTMB::betabinomial(link = "logit"),
  control = glmmTMB::glmmTMBControl(
    optimizer = optim, 
    optArgs = list(method = "BFGS")
  )
)
#> Warning in eval(family$initialize): non-integer #successes in a binomial
#> glm!
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model
#> convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')

mbepr_model
#> Formula:          mbepr ~ spei_12m + gdp_per_capita + (1 | year)
#> Data: data
#>      AIC      BIC   logLik df.resid 
#>       NA       NA       NA    54665 
#> Random-effects (co)variances:
#> 
#> Conditional model:
#>  Groups Name        Std.Dev.
#>  year   (Intercept) 1.512832
#> 
#> Number of obs: 54670 / Conditional model: year, 10
#> 
#> Dispersion parameter for betabinomial family (): 5.2e+177 
#> 
#> Fixed Effects:
#> 
#> Conditional model:
#>      (Intercept)          spei_12m    gdp_per_capita  
#> -0.5702246509759  -5.7450764201136   0.0000008241382
Code
overdispersion_test <- 
  mbepr_model |>
  stats::residuals(type = "pearson") %>%
  `^`(2) %>%
  `/`(mbepr_model$df.residual) |>
  sum()
cli::cli_inform("Overdispersion statistic: {overdispersion_test}\n")
#> Overdispersion statistic: 0

if (overdispersion_test > 1.5) {
  cli::cli_inform("There might still be overdispersion in the model.\n")
}
Code
dplyr::tibble(
  fitted = stats::fitted(mbepr_model),
  residuals = stats::residuals(mbepr_model, type = "pearson")
) |>
  ggplot2::ggplot(ggplot2::aes(x = fitted, y = residuals)) +
  ggplot2::geom_point(alpha = 0.5) +
  ggplot2::geom_hline(
    yintercept = 0, 
    color = get_brand_color("red")
  ) +
  ggplot2::labs(
    x = "Fitted values",
    y = "Pearson residuals",
    title = "Residuals vs. Fitted"
  )

Code
data |>
  ggplot2::ggplot(ggplot2::aes(x = spei_12m, y = mbepr)) +
  ggplot2::geom_point(alpha = 0.5) +
  ggplot2::geom_smooth(
    method = "glm", 
    method.args = list(family = "binomial"), 
    se = TRUE, 
    color = get_brand_color("red")
  ) +
  ggplot2::labs(
    x = "SPEI (12 months)",
    y = "Probability of MBEPR",
    title = "Relationship between SPEI_12m and MBEPR"
  )
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning in eval(family$initialize): non-integer #successes in a binomial
#> glm!

Code
data |>
  tidyr::drop_na(gdp_per_capita, mbepr) |>
  ggplot2::ggplot(ggplot2::aes(x = gdp_per_capita, y = mbepr)) +
  ggplot2::geom_point(alpha = 0.5) +
  ggplot2::geom_smooth(
    method = "glm", 
    method.args = list(family = "binomial"), 
    se = TRUE, 
    color = get_brand_color("red")
  ) +
  ggplot2::labs(
    x = "GDP per capita",
    y = "Very Short Stature for Age (MBEPR)",
    title = "Relationship between GDP per capita and MBEPR"
  )
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning in eval(family$initialize): non-integer #successes in a binomial
#> glm!

BEIPR – Short Stature for Age (Baixa Estatura Para Idade)

beipr_model <- glmmTMB::glmmTMB(
  beipr ~ spei_12m * gdp_per_capita + (1 | year),
  data = data,
  family = glmmTMB::betabinomial(link = "logit")
)
#> Warning in eval(family$initialize): non-integer #successes in a binomial
#> glm!
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model
#> convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model
#> convergence problem; false convergence (8). See vignette('troubleshooting'),
#> help('diagnose')

beipr_model
#> Formula:          beipr ~ spei_12m * gdp_per_capita + (1 | year)
#> Data: data
#>      AIC      BIC   logLik df.resid 
#>       NA       NA       NA    54664 
#> Random-effects (co)variances:
#> 
#> Conditional model:
#>  Groups Name        Std.Dev.    
#>  year   (Intercept) 0.0006285806
#> 
#> Number of obs: 54670 / Conditional model: year, 10
#> 
#> Dispersion parameter for betabinomial family (): 2.47e+08 
#> 
#> Fixed Effects:
#> 
#> Conditional model:
#>             (Intercept)                 spei_12m           gdp_per_capita  
#>        -2.5915189700578         -0.0470766191722         -0.0000050110971  
#> spei_12m:gdp_per_capita  
#>        -0.0000004649216
summary(beipr_model)
#>  Family: betabinomial  ( logit )
#> Formula:          beipr ~ spei_12m * gdp_per_capita + (1 | year)
#> Data: data
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>       NA       NA       NA       NA    54664 
#> 
#> Random effects:
#> 
#> Conditional model:
#>  Groups Name        Variance        Std.Dev.    
#>  year   (Intercept) 0.0000003951135 0.0006285806
#> Number of obs: 54670, groups:  year, 10
#> 
#> Dispersion parameter for betabinomial family (): 2.47e+08 
#> 
#> Conditional model:
#>                                 Estimate       Std. Error    z value
#> (Intercept)             -2.5915189700578  0.0245448649134 -105.58294
#> spei_12m                -0.0470766191722  0.0881768574492   -0.53389
#> gdp_per_capita          -0.0000050110971  0.0000008757264   -5.72222
#> spei_12m:gdp_per_capita -0.0000004649216  0.0000053852446   -0.08633
#>                               Pr(>|z|)    
#> (Intercept)                 < 2.22e-16 ***
#> spei_12m                       0.59342    
#> gdp_per_capita          0.000000010514 ***
#> spei_12m:gdp_per_capita        0.93120    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MAPER – Severe Thinness for Height (Magreza Acentuada Para a Estatura)

maper_model <- glm(
  maper ~ spei_12m * gdp_per_capita,
  data = data,
  family = quasibinomial(link = "logit")
)

maper_model
#> 
#> Call:  glm(formula = maper ~ spei_12m * gdp_per_capita, family = quasibinomial(link = "logit"), 
#>     data = data)
#> 
#> Coefficients:
#>             (Intercept)                 spei_12m           gdp_per_capita  
#>         -3.536476647077          -0.057882158222          -0.000011753764  
#> spei_12m:gdp_per_capita  
#>         -0.000005560497  
#> 
#> Degrees of Freedom: 54669 Total (i.e. Null);  54666 Residual
#>   (10331 observations deleted due to missingness)
#> Null Deviance:       1046.018 
#> Residual Deviance: 1008.234  AIC: NA
summary(maper_model)
#> 
#> Call:
#> glm(formula = maper ~ spei_12m * gdp_per_capita, family = quasibinomial(link = "logit"), 
#>     data = data)
#> 
#> Coefficients:
#>                                 Estimate       Std. Error    t value
#> (Intercept)             -3.5364766470774  0.0106686041904 -331.48447
#> spei_12m                -0.0578821582221  0.0211670993397   -2.73453
#> gdp_per_capita          -0.0000117537643  0.0000004997133  -23.52102
#> spei_12m:gdp_per_capita -0.0000055604975  0.0000010585039   -5.25317
#>                              Pr(>|t|)    
#> (Intercept)                < 2.22e-16 ***
#> spei_12m                    0.0062489 ** 
#> gdp_per_capita             < 2.22e-16 ***
#> spei_12m:gdp_per_capita 0.00000015006 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for quasibinomial family taken to be 0.0366800415)
#> 
#>     Null deviance: 1046.0178  on 54669  degrees of freedom
#> Residual deviance: 1008.2344  on 54666  degrees of freedom
#>   (10331 observations deleted due to missingness)
#> AIC: NA
#> 
#> Number of Fisher Scoring iterations: 7

MPEPR – Thinness for Height (Magreza Por Estatura)

mpepr_model <- glm(
  mpepr ~ spei_12m * gdp_per_capita,
  data = data,
  family = quasibinomial(link = "logit")
)

mpepr_model
#> 
#> Call:  glm(formula = mpepr ~ spei_12m * gdp_per_capita, family = quasibinomial(link = "logit"), 
#>     data = data)
#> 
#> Coefficients:
#>             (Intercept)                 spei_12m           gdp_per_capita  
#>         -3.453653181925          -0.129070860970          -0.000007061586  
#> spei_12m:gdp_per_capita  
#>         -0.000002473331  
#> 
#> Degrees of Freedom: 54669 Total (i.e. Null);  54666 Residual
#>   (10331 observations deleted due to missingness)
#> Null Deviance:       687.2864 
#> Residual Deviance: 661.839   AIC: NA
summary(mpepr_model)
#> 
#> Call:
#> glm(formula = mpepr ~ spei_12m * gdp_per_capita, family = quasibinomial(link = "logit"), 
#>     data = data)
#> 
#> Coefficients:
#>                                 Estimate       Std. Error    t value
#> (Intercept)             -3.4536531819248  0.0055451302278 -622.82634
#> spei_12m                -0.1290708609697  0.0110714003091  -11.65804
#> gdp_per_capita          -0.0000070615864  0.0000002407524  -29.33132
#> spei_12m:gdp_per_capita -0.0000024733309  0.0000005202559   -4.75407
#>                             Pr(>|t|)    
#> (Intercept)               < 2.22e-16 ***
#> spei_12m                  < 2.22e-16 ***
#> gdp_per_capita            < 2.22e-16 ***
#> spei_12m:gdp_per_capita 0.0000019987 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for quasibinomial family taken to be 0.01240486978)
#> 
#>     Null deviance: 687.28635  on 54669  degrees of freedom
#> Residual deviance: 661.83899  on 54666  degrees of freedom
#>   (10331 observations deleted due to missingness)
#> AIC: NA
#> 
#> Number of Fisher Scoring iterations: 7

Generalized Additive Model Mixed Effects (GAMM)

MBEPR – Very Short Stature for Age (Muito Baixa Estatura Para a Idade)

Code
data <-
  data |>
  dplyr::mutate(
    mbepr = pmax(0.0001, pmin(0.9999, mbepr))
  )
Code
mbepr_model_gamm <- mgcv::gam(
  mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re"), 
  data = data,
  family = mgcv::betar(link = "logit")
)

mbepr_model_gamm
#> 
#> Family: Beta regression(24.531) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Estimated degrees of freedom:
#> 8.89 0.00  total = 10.89 
#> 
#> REML score: -106831.1863
Code
summary(mbepr_model_gamm)
#> 
#> Family: Beta regression(24.531) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -2.7027669847549  0.0063418020582 -426.18280 < 2.22e-16 ***
#> gdp_per_capita -0.0000086801261  0.0000001794867  -48.36084 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df    Chi.sq    p-value    
#> s(spei_12m) 8.889552237951 8.996491 928.90218 < 2.22e-16 ***
#> s(year)     0.000003234013 1.000000   0.00005 0.00019065 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0243   Deviance explained = 5.05%
#> -REML = -1.0683e+05  Scale est. = 1         n = 54670
Code
mbepr_model_gamm |> plot(pages = 1)

Code
mbepr_model_gamm |> coef_mean()
#> [1] -0.2619787413
Code
mbepr_model_gamm |> coef_spei12m_mean()
#> [1] -0.04899657958
Code
data |>
plot_gamm(
  model = mbepr_model_gamm,
  title = paste0(
    "Very Short Stature for Age versus (MBEPR) versus \n", 
    "Standardised Precipitation Evapotranspiration Index"
  ),
  y_label = "Predicted MBEPR"
)

BEIPR – Short Stature for Age (Baixa Estatura Para Idade)

Code
data <-
  data |>
  dplyr::mutate(
    beipr = pmax(0.0001, pmin(0.9999, beipr))
  )
Code
beipr_model_gamm <- mgcv::gam(
  beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re"),
  data = data,
  family = mgcv::betar(link = "logit")
)

beipr_model_gamm
#> 
#> Family: Beta regression(35.93) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Estimated degrees of freedom:
#> 8.69 0.00  total = 10.69 
#> 
#> REML score: -105270.0952
Code
summary(beipr_model_gamm)
#> 
#> Family: Beta regression(35.93) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -2.5588107330787  0.0083680633557 -305.78291 < 2.22e-16 ***
#> gdp_per_capita -0.0000065235994  0.0000001544059  -42.24967 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df    Chi.sq    p-value    
#> s(spei_12m) 8.68890721298 8.972917 679.26856 < 2.22e-16 ***
#> s(year)     0.00001312662 1.000000   0.00081 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0283   Deviance explained = 4.13%
#> -REML = -1.0527e+05  Scale est. = 1         n = 54670
Code
beipr_model_gamm |> plot(pages = 1)

Code
beipr_model_gamm |> coef_mean()
#> [1] -0.2643808837
Code
beipr_model_gamm |> coef_spei12m_mean()
#> [1] -0.06819482812
Code
data |>
plot_gamm(
  model = beipr_model_gamm,
  title = paste0(
    "Short Stature for Age (BEIPR) versues \n",
    "Standardised Precipitation Evapotranspiration Index"
  ),
  x_label = "SPEI 12m",
  y_label = "Predicted BEIPR"
)

MAPER – Severe Thinness for Height (Magreza Acentuada Para a Estatura)

Code
data <-
  data |>
  dplyr::mutate(
    maper = pmax(0.0001, pmin(0.9999, maper))
  )
Code
maper_model_gamm <- mgcv::gam(
  maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re"),
  data = data,
  family = mgcv::betar(link = "logit")
)

maper_model_gamm
#> 
#> Family: Beta regression(37.008) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Estimated degrees of freedom:
#> 8.86 0.00  total = 10.86 
#> 
#> REML score: -151474.9192
Code
summary(maper_model_gamm)
#> 
#> Family: Beta regression(37.008) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -3.5187320870973  0.0065372134986 -538.26177 < 2.22e-16 ***
#> gdp_per_capita -0.0000110852034  0.0000001906241  -58.15216 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df     Chi.sq    p-value    
#> s(spei_12m) 8.864908270421 8.994767 1385.13687 < 2.22e-16 ***
#> s(year)     0.000002548998 1.000000    0.00013 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0272   Deviance explained =  6.9%
#> -REML = -1.5147e+05  Scale est. = 1         n = 54670
Code
maper_model_gamm |> plot(pages = 1)

Code
maper_model_gamm |> coef_mean()
#> [1] -0.3810039149
Code
maper_model_gamm |> coef_spei12m_mean()
#> [1] -0.1170337536
Code
data |>
plot_gamm(
  model = maper_model_gamm,
  title = paste0(
    "Severe Thinness for Height (MAPER) versus \n",
    "Standardised Precipitation Evapotranspiration Index"
  ),
  y_label = "Predicted MAPER"
)

MPEPR – Thinness for Height (Magreza Por Estatura)

Code
data <-
  data |>
  dplyr::mutate(
    mpepr = pmax(0.0001, pmin(0.9999, mpepr))
  )
Code
mpepr_model_gamm <- mgcv::gam(
  mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re"),
  data = data,
  family = mgcv::betar(link = "logit")
)

mpepr_model_gamm
#> 
#> Family: Beta regression(46.805) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Estimated degrees of freedom:
#> 8.88 0.00  total = 10.88 
#> 
#> REML score: -143159.3615
Code
summary(mpepr_model_gamm)
#> 
#> Family: Beta regression(46.805) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error   z value   Pr(>|z|)    
#> (Intercept)    -3.3850047188458  0.0062439242831 -542.1278 < 2.22e-16 ***
#> gdp_per_capita -0.0000094624353  0.0000001793694  -52.7539 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df   Chi.sq p-value    
#> s(spei_12m) 8.88155042288 8.995962 2028.904 < 2e-16 ***
#> s(year)     0.00000319169 1.000000    0.000 0.46573    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0337   Deviance explained = 7.53%
#> -REML = -1.4316e+05  Scale est. = 1         n = 54670
Code
mpepr_model_gamm |> plot(pages = 1)

Code
mpepr_model_gamm |> coef_mean()
#> [1] -0.4039970943
Code
mpepr_model_gamm |> coef_spei12m_mean()
#> [1] -0.1625501059
Code
data |>
plot_gamm(
  model = mpepr_model_gamm,
  title = paste0(
    "Thinness for Height (MPEPR) versus \n",
    "Standardised Precipitation Evapotranspiration Index"
  ),
  y_label = "Predicted MPEPR"
)

Model Check

Code
mgcv::gam.check(mbepr_model_gamm)

#> 
#> Method: REML   Optimizer: outer newton
#> full convergence after 5 iterations.
#> Gradient range [0.00002093215649,0.09068848415]
#> (score -106831.1863 & scale 1).
#> eigenvalue range [-0.00002093202182,24100.629].
#> Model rank =  12 / 12 
#> 
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#> 
#>                     k'        edf k-index p-value  
#> s(spei_12m) 9.00000000 8.88955224    0.99   0.610  
#> s(year)     1.00000000 0.00000323    0.97   0.045 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
mgcv::gam.check(beipr_model_gamm)

#> 
#> Method: REML   Optimizer: outer newton
#> full convergence after 5 iterations.
#> Gradient range [-0.000215252209,0.0003964991776]
#> (score -105270.0952 & scale 1).
#> eigenvalue range [-0.0003964890913,25913.25923].
#> Model rank =  12 / 12 
#> 
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#> 
#>                    k'       edf k-index p-value
#> s(spei_12m) 9.0000000 8.6889072    1.01    0.85
#> s(year)     1.0000000 0.0000131    0.98    0.17
Code
mgcv::gam.check(maper_model_gamm)

#> 
#> Method: REML   Optimizer: outer newton
#> full convergence after 4 iterations.
#> Gradient range [0.00006354959436,0.01155840506]
#> (score -151474.9192 & scale 1).
#> eigenvalue range [-0.00006354927352,21016.01011].
#> Model rank =  12 / 12 
#> 
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#> 
#>                     k'        edf k-index p-value
#> s(spei_12m) 9.00000000 8.86490827    1.01    0.71
#> s(year)     1.00000000 0.00000255    1.02    0.86
Code
mgcv::gam.check(mpepr_model_gamm)

#> 
#> Method: REML   Optimizer: outer newton
#> full convergence after 4 iterations.
#> Gradient range [-0.000000700510767,0.01755277868]
#> (score -143159.3615 & scale 1).
#> Hessian positive definite, eigenvalue range [0.0000007005112393,23265.15506].
#> Model rank =  12 / 12 
#> 
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#> 
#>                     k'        edf k-index p-value
#> s(spei_12m) 9.00000000 8.88155042    1.01    0.89
#> s(year)     1.00000000 0.00000319    0.98    0.25

Analysis by MISFS

By S(SPEI 12m) + GDP Per Capita + S(Year)

MBEPR – Very Short Stature for Age (Muito Baixa Estatura Para a Idade)

Code
mbepr_model_gamm_misfs_1 <- 
  data |> 
  model_gamm_misfs(
    model = mbepr_model_gamm,
    type = 1
  )
mbepr_model_gamm_misfs_1 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(24.413) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value         Pr(>|z|)
#> (Intercept)    -2.7086722214888  0.0251079118483 -107.88122       < 2.22e-16
#> gdp_per_capita -0.0000027789004  0.0000004522482   -6.14464 0.00000000080148
#>                   
#> (Intercept)    ***
#> gdp_per_capita ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df   Chi.sq     p-value    
#> s(spei_12m) 6.907630259034 8.023555 80.29324  < 2.22e-16 ***
#> s(year)     0.000009175927 1.000000  0.00014 0.000077048 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00638   Deviance explained = 1.88%
#> -REML = -12621  Scale est. = 1         n = 6721
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(22.125) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                       Estimate      Std. Error   z value   Pr(>|z|)    
#> (Intercept)    -3.035549328216  0.009437220353 -321.6571 < 2.22e-16 ***
#> gdp_per_capita -0.000001868094  0.000000216192   -8.6409 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df    Chi.sq    p-value    
#> s(spei_12m) 8.185709219345 8.826234 114.10441 < 2.22e-16 ***
#> s(year)     0.000002647121 1.000000   0.00002  0.0028869 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.00115   Deviance explained = 0.636%
#> -REML = -59428  Scale est. = 1         n = 27944
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(44.041) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value Pr(>|z|)    
#> (Intercept)    -2.6498919545521  0.0144212926778 -183.74857   <2e-16 ***
#> gdp_per_capita -0.0000007112901  0.0000004455574   -1.59641   0.1104    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf  Ref.df    Chi.sq    p-value    
#> s(spei_12m) 8.79678710055 8.98819 205.55764 < 2.22e-16 ***
#> s(year)     0.00001311438 1.00000   0.00131 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0121   Deviance explained = 1.18%
#> -REML = -35016  Scale est. = 1         n = 17659
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(39.046) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error   z value   Pr(>|z|)    
#> (Intercept)    -2.2130292815722  0.0398572683822 -55.52386 < 2.22e-16 ***
#> gdp_per_capita -0.0000088877279  0.0000009970597  -8.91394 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df   Chi.sq    p-value    
#> s(spei_12m) 1.00108716051 1.002174 11.00261 0.00091716 ***
#> s(year)     0.00002191564 1.000000  0.00023 0.00115144 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.034   Deviance explained =    4%
#> -REML = -4102.3  Scale est. = 1         n = 2346
mbepr_model_gamm_misfs_1 |> coef_means_misfs()
mbepr_model_gamm_misfs_1 |> coef_spei12m_means_misfs()
Code
data |>
  plot_gamm_misfs(
    gamm_models = mbepr_model_gamm_misfs_1,
    title = "Very Short Stature for Age (MBEPR) by MISFS",
    y_label = "Predicted MBEPR",
    type = 1
  )

BEIPR – Short Stature for Age (Baixa Estatura Para Idade)

Code
beipr_model_gamm_misfs_1 <- 
  data |> 
  model_gamm_misfs(
    model = beipr_model_gamm,
    type = 1
  )
beipr_model_gamm_misfs_1 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(33.872) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value    Pr(>|z|)    
#> (Intercept)    -2.6183178600095  0.0242319981512 -108.05208  < 2.22e-16 ***
#> gdp_per_capita -0.0000016539806  0.0000003813088   -4.33764 0.000014402 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df    Chi.sq    p-value    
#> s(spei_12m) 7.84875841922 8.666477 161.36317 < 2.22e-16 ***
#> s(year)     0.00001262858 1.000000   0.00115 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0151   Deviance explained = 2.68%
#> -REML = -12689  Scale est. = 1         n = 6721
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(32.684) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value         Pr(>|z|)
#> (Intercept)    -2.8323627140074  0.0086247719049 -328.39856       < 2.22e-16
#> gdp_per_capita -0.0000011254321  0.0000001797438   -6.26131 0.00000000038175
#>                   
#> (Intercept)    ***
#> gdp_per_capita ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df   Chi.sq      p-value    
#> s(spei_12m) 4.169516002121 5.191392 55.19505   < 2.22e-16 ***
#> s(year)     0.000003964537 1.000000  0.00010 0.0000008686 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.00232   Deviance explained = 0.34%
#> -REML = -56179  Scale est. = 1         n = 27944
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(82.194) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -2.5322023592576  0.0091352621940 -277.18989 < 2.22e-16 ***
#> gdp_per_capita -0.0000012765615  0.0000003306047   -3.86129 0.00011279 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df    Chi.sq    p-value    
#> s(spei_12m) 8.662526874889 8.968023 121.74172 < 2.22e-16 ***
#> s(year)     0.000009311943 1.000000   0.00029 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0106   Deviance explained = 0.802%
#> -REML = -38558  Scale est. = 1         n = 17659
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(59.763) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error   z value   Pr(>|z|)    
#> (Intercept)    -1.8656879394534  0.0394343023213 -47.31130 < 2.22e-16 ***
#> gdp_per_capita -0.0000072831416  0.0000007320574  -9.94887 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df   Chi.sq   p-value   
#> s(spei_12m) 1.00047904420 1.000958 10.21630 0.0013953 **
#> s(year)     0.00004385054 1.000000  0.00026 0.0158747 * 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0508   Deviance explained = 4.65%
#> -REML =  -4141  Scale est. = 1         n = 2346
beipr_model_gamm_misfs_1 |> coef_means_misfs()
beipr_model_gamm_misfs_1 |> coef_spei12m_means_misfs()
Code
data |>
  plot_gamm_misfs(
    gamm_models = beipr_model_gamm_misfs_1,
    title = "Short Stature for Age (BEIPR) by MISFS",
    y_label = "Predicted BEIPR",
    type = 1
  )

MAPER – Severe Thinness for Height (Magreza Acentuada Para a Estatura)

Code
maper_model_gamm_misfs_1 <- 
  data |> 
  model_gamm_misfs(
    model = maper_model_gamm,
    type = 1
  )
maper_model_gamm_misfs_1 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(33.407) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error   z value        Pr(>|z|)
#> (Intercept)    -3.4874339100552  0.0195070801919 -178.7779      < 2.22e-16
#> gdp_per_capita -0.0000029806451  0.0000004977948   -5.9877 0.0000000021283
#>                   
#> (Intercept)    ***
#> gdp_per_capita ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df   Chi.sq     p-value    
#> s(spei_12m) 5.460870605461 6.653486 32.32437 0.000030034 ***
#> s(year)     0.000002606773 1.000000  0.00002   0.0051981 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0035   Deviance explained = 1.14%
#> -REML = -17373  Scale est. = 1         n = 6721
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(33.321) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -3.8890807401487  0.0098061169347 -396.59743 < 2.22e-16 ***
#> gdp_per_capita -0.0000033276989  0.0000002339366  -14.22479 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df    Chi.sq  p-value    
#> s(spei_12m) 8.13617255244 8.805671 246.73101  < 2e-16 ***
#> s(year)     0.00000214972 1.000000   0.00001 0.019722 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00601   Deviance explained =  1.4%
#> -REML = -86143  Scale est. = 1         n = 27944
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(62.482) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                       Estimate      Std. Error    z value Pr(>|z|)    
#> (Intercept)    -3.436883594232  0.010291548864 -333.95203  < 2e-16 ***
#> gdp_per_capita -0.000000542129  0.000000507238   -1.06879  0.28517    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df    Chi.sq    p-value    
#> s(spei_12m) 6.936506246539 8.043066 120.57074 < 2.22e-16 ***
#> s(year)     0.000003473506 1.000000   0.00047 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00537   Deviance explained = 0.755%
#> -REML = -45685  Scale est. = 1         n = 17659
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(96.599) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error   z value     Pr(>|z|)    
#> (Intercept)    -3.4933442183346  0.0384152852640 -90.93631   < 2.22e-16 ***
#> gdp_per_capita -0.0000046079335  0.0000009507086  -4.84684 0.0000012544 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf  Ref.df  Chi.sq      p-value    
#> s(spei_12m) 1.23786573662 1.44072 0.12676      0.91803    
#> s(year)     0.00001716901 1.00000 0.00042 0.0000013159 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00818   Deviance explained = 1.16%
#> -REML = -6548.9  Scale est. = 1         n = 2346
maper_model_gamm_misfs_1 |> coef_means_misfs()
maper_model_gamm_misfs_1 |> coef_spei12m_means_misfs()
Code
data |>
  plot_gamm_misfs(
    gamm_models = maper_model_gamm_misfs_1,
    title = "Severe Thinness for Height (MAPER) by MISFS",
    y_label = "Predicted MAPER",
    type = 1
  )

MPEPR – Thinness for Height (Magreza Por Estatura)

Code
mpepr_model_gamm_misfs_1 <- 
  data |> 
  model_gamm_misfs(
    model = mpepr_model_gamm,
    type = 1
  )
mpepr_model_gamm_misfs_1 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(44.225) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value     Pr(>|z|)    
#> (Intercept)    -3.3687440829323  0.0246335647481 -136.75423   < 2.22e-16 ***
#> gdp_per_capita -0.0000021283608  0.0000004427547   -4.80709 0.0000015314 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df   Chi.sq p-value    
#> s(spei_12m) 5.970719189065 7.173862 65.21956  <2e-16 ***
#> s(year)     0.000009043833 1.000000  0.00002  0.1332    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00253   Deviance explained = 1.41%
#> -REML = -16700  Scale est. = 1         n = 6721
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(37.503) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -3.6806629571423  0.0094521649778 -389.39893 < 2.22e-16 ***
#> gdp_per_capita -0.0000025259286  0.0000002208273  -11.43848 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df    Chi.sq     p-value    
#> s(spei_12m) 8.579675779618 8.950963 521.69511  < 2.22e-16 ***
#> s(year)     0.000002501037 1.000000   0.00005 0.000010739 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00796   Deviance explained = 2.15%
#> -REML = -77547  Scale est. = 1         n = 27944
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(109.854) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value Pr(>|z|)    
#> (Intercept)    -3.3387629360220  0.0138461644147 -241.13269  < 2e-16 ***
#> gdp_per_capita -0.0000011686712  0.0000003990934   -2.92831 0.003408 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                       edf   Ref.df   Chi.sq    p-value    
#> s(spei_12m) 7.38504143582 8.380703 97.23446 < 2.22e-16 ***
#> s(year)     0.00001612837 1.000000  0.00249 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00699   Deviance explained = 0.647%
#> -REML = -48207  Scale est. = 1         n = 17659
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(128.082) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(spei_12m) + gdp_per_capita + s(year, bs = "re")
#> 
#> Parametric coefficients:
#>                        Estimate       Std. Error    z value   Pr(>|z|)    
#> (Intercept)    -3.3036765685165  0.0246968794470 -133.76899 < 2.22e-16 ***
#> gdp_per_capita -0.0000019221267  0.0000006966528   -2.75909  0.0057963 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>                        edf   Ref.df  Chi.sq    p-value    
#> s(spei_12m) 1.002424084109 1.004844 0.08108 0.77843854    
#> s(year)     0.000009718337 1.000000 0.00013 0.00027477 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.00351   Deviance explained = 0.34%
#> -REML = -6506.5  Scale est. = 1         n = 2346
mpepr_model_gamm_misfs_1 |> coef_means_misfs()
mpepr_model_gamm_misfs_1 |> coef_spei12m_means_misfs()
Code
data |>
  plot_gamm_misfs(
    gamm_models = mpepr_model_gamm_misfs_1,
    title = "Thinness for Height (MPEPR) by MISFS",
    y_label = "Predicted MPEPR",
    type = 1
  )

Only by S(Year)

MBEPR – Very Short Stature for Age (Muito Baixa Estatura Para a Idade)

Code
mbepr_model_gamm_misfs_2 <- 
  data |> 
  model_gamm_misfs(
    model = mbepr_model_gamm,
    type = 2
  )
mbepr_model_gamm_misfs_2 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(21.813) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -2.760548281  0.008408097 -328.3202 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 7.636547 8.537803 109.6488 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.0117   Deviance explained = 1.55%
#> -REML = -14802  Scale est. = 1         n = 8001
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(20.648) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.068052812  0.004517026 -679.2197 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf  Ref.df   Chi.sq    p-value    
#> s(year) 7.967059 8.72326 184.8139 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.00213   Deviance explained = 0.595%
#> -REML = -70253  Scale est. = 1         n = 33201
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(39.594) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(year)
#> 
#> Parametric coefficients:
#>                Estimate  Std. Error   z value   Pr(>|z|)    
#> (Intercept) -2.61385934  0.00398162 -656.4814 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>             edf   Ref.df   Chi.sq    p-value    
#> s(year) 8.13964 8.803401 449.4811 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0303   Deviance explained = 2.12%
#> -REML = -40355  Scale est. = 1         n = 21009
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(35.132) 
#> Link function: logit 
#> 
#> Formula:
#> mbepr ~ s(year)
#> 
#> Parametric coefficients:
#>                Estimate  Std. Error   z value   Pr(>|z|)    
#> (Intercept) -2.26649295  0.01021249 -221.9333 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf  Ref.df   Chi.sq    p-value    
#> s(year) 6.616444 7.75189 125.4656 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0429   Deviance explained = 4.53%
#> -REML = -4698.8  Scale est. = 1         n = 2790
mbepr_model_gamm_misfs_2 |> coef_means_misfs()
Code
summarise_model_gamm_misfs_2(mbepr_model_gamm_misfs_2)
Code
data |>
  plot_gamm_misfs(
    gamm_models = mbepr_model_gamm_misfs_2,
    title = "Very Short Stature for Age (MBEPR) by MISFS",
    x_label = "Years",
    y_label = "Predicted MBEPR",
    type = 2
  )

BEIPR – Short Stature for Age (Baixa Estatura Para Idade)

Code
beipr_model_gamm_misfs_2 <- 
  data |> 
  model_gamm_misfs(
    model = beipr_model_gamm,
    type = 2
  )
beipr_model_gamm_misfs_2 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(29.769) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error  z value   Pr(>|z|)    
#> (Intercept) -2.648751917  0.007303033 -362.692 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df  Chi.sq    p-value    
#> s(year) 6.529559 7.673822 147.148 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.0162   Deviance explained = 1.92%
#> -REML = -14767  Scale est. = 1         n = 8001
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(29.986) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -2.847493774  0.003799504 -749.4383 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf  Ref.df   Chi.sq    p-value    
#> s(year) 8.066288 8.77101 177.0859 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.00791   Deviance explained = 0.587%
#> -REML = -65744  Scale est. = 1         n = 33201
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(74.565) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(year)
#> 
#> Parametric coefficients:
#>                Estimate  Std. Error   z value   Pr(>|z|)    
#> (Intercept) -2.51098479  0.00290845 -863.3411 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 7.426274 8.398473 562.5951 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0382   Deviance explained = 2.59%
#> -REML = -44659  Scale est. = 1         n = 21009
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(52.489) 
#> Link function: logit 
#> 
#> Formula:
#> beipr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -1.924461480  0.007582017 -253.8192 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 3.601066 4.457232 66.45076 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.032   Deviance explained = 2.47%
#> -REML = -4742.4  Scale est. = 1         n = 2790
beipr_model_gamm_misfs_2 |> coef_means_misfs()
Code
summarise_model_gamm_misfs_2(beipr_model_gamm_misfs_2)
Code
data |>
  plot_gamm_misfs(
    gamm_models = beipr_model_gamm_misfs_2,
    title = "Short Stature for Age (BEIPR) by MISFS",
    y_label = "Predicted BEIPR",
    type = 2
  )

MAPER – Severe Thinness for Height (Magreza Acentuada Para a Estatura)

Code
maper_model_gamm_misfs_2 <- 
  data |> 
  model_gamm_misfs(
    model = maper_model_gamm,
    type = 2
  )
maper_model_gamm_misfs_2 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(30.163) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.528075401  0.009245944 -381.5809 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 6.073143 7.239167 43.63848 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.00192   Deviance explained = 0.64%
#> -REML = -20507  Scale est. = 1         n = 8001
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(30.916) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.939262865  0.004881495 -806.9787 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 8.213052 8.834218 170.5708 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.000345   Deviance explained = 0.525%
#> -REML = -1.019e+05  Scale est. = 1         n = 33201
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(57.111) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.411791854  0.004535189 -752.2932 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 8.341086 8.881764 269.4013 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0154   Deviance explained = 1.38%
#> -REML = -53406  Scale est. = 1         n = 21009
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(90.106) 
#> Link function: logit 
#> 
#> Formula:
#> maper ~ s(year)
#> 
#> Parametric coefficients:
#>                Estimate  Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.51658361  0.01083901 -324.4378 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 3.361261 4.168493 55.80129 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0261   Deviance explained = 2.21%
#> -REML =  -7679  Scale est. = 1         n = 2790
maper_model_gamm_misfs_2 |> coef_means_misfs()
Code
summarise_model_gamm_misfs_2(maper_model_gamm_misfs_2)
Code
data |>
  plot_gamm_misfs(
    gamm_models = maper_model_gamm_misfs_2,
    title = "Severe Thinness for Height (MAPER) by MISFS",
    y_label = "Predicted MAPER",
    type = 2
  )

MPEPR – Thinness for Height (Magreza Por Estatura)

Code
mpepr_model_gamm_misfs_2 <- 
  data |> 
  model_gamm_misfs(
    model = mpepr_model_gamm,
    type = 2
  )
mpepr_model_gamm_misfs_2 |> summarise_model_gamm_misfs()
#> 
#> ── Model A ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(39.693) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.412541413  0.008316807 -410.3187 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 5.513303 6.660799 92.81762 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.0187   Deviance explained = 1.23%
#> -REML = -19723  Scale est. = 1         n = 8001
#> 
#> ── Model B ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(34.92) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.725866332  0.004572887 -814.7733 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df  Chi.sq    p-value    
#> s(year) 8.256788 8.851361 335.879 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  -0.00327   Deviance explained = 0.97%
#> -REML = -91740  Scale est. = 1         n = 33201
#> 
#> ── Model C ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(98.152) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(year)
#> 
#> Parametric coefficients:
#>                 Estimate   Std. Error   z value   Pr(>|z|)    
#> (Intercept) -3.335865028  0.003554644 -938.4527 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq    p-value    
#> s(year) 5.922402 7.086925 151.7835 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0133   Deviance explained = 0.766%
#> -REML = -56263  Scale est. = 1         n = 21009
#> 
#> ── Model D ──────────────────────────────────────────────────────────────────
#> 
#> Family: Beta regression(123.387) 
#> Link function: logit 
#> 
#> Formula:
#> mpepr ~ s(year)
#> 
#> Parametric coefficients:
#>                Estimate  Std. Error  z value   Pr(>|z|)    
#> (Intercept) -3.31792430  0.00876765 -378.428 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf   Ref.df   Chi.sq     p-value    
#> s(year) 2.625175 3.266502 22.89967 0.000070824 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.0108   Deviance explained = 0.925%
#> -REML = -7695.9  Scale est. = 1         n = 2790
mpepr_model_gamm_misfs_2 |> coef_means_misfs()
Code
summarise_model_gamm_misfs_2(mpepr_model_gamm_misfs_2)
Code
data |>
  plot_gamm_misfs(
    gamm_models = mpepr_model_gamm_misfs_2,
    title = "Thinness for Height (MPEPR) by MISFS",
    y_label = "Predicted MPEPR",
    type = 2
  )

Acknowledgments

This analysis is part of the Sustentarea Research and Extension Group’s project: Global syndemic: the impact of anthropogenic climate change on the health and nutrition of children under five years old attended by Brazil’s public health system (SUS).

This research was supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brazil (CNPq).

References

Wickham, H., Çetinkaya-Rundel, M., & Grolemund, G. (2023). R for data science: Import, tidy, transform, visualize, and model data (2nd ed.). O’Reilly Media. https://r4ds.hadley.nz